//
//  RandomTree.cpp
//  RandomForest-AndreaBalboni
//
//  Created by Andrea Balboni on 05/04/12.
//  Copyright (c) 2012 Università degli Studi di Modena e Reggio Emilia. All rights reserved.
//
#include <iostream>
#include <algorithm>
#include <vector>
#include "RandomTree.h"
#include <opencv2/opencv.hpp>


#define RANDOMIZE_GINI_INDEX

/**
 * Comparator functors
 */

/**
 * Funzione template usata come comparatore per l'ordinamento del vettore delle feature per il gini index 
 */
template <typename _T> bool featureComparator(std::pair<_T, unsigned int> el1, std::pair<_T, unsigned int> el2) {
	return el1.first < el2.first;
}

bool classProbComparator(std::pair<int,float> c1, std::pair<int, float> c2) {
    return c1.second<c2.second;
}

template <typename _T>
bool distComparator(std::pair<float,_T> c1, std::pair<float, _T> c2) {
    return c1.first<c2.first;
}

/* 
 * Utility functions
 */

template <typename _T>
 void RandomTree<_T>::getNodeClassMap(std::vector<unsigned int> &indexes, TrainingData<_T> &DataSet, std::map<int, float> &classMap){
    std::vector<unsigned int>::iterator it = indexes.begin();
    std::map<int,float>::iterator mapIter;
    for(;it!=indexes.end();++it){
        int classe=DataSet.getRecord(*it).getClass();
        mapIter=classMap.find(classe);
        if(mapIter!=classMap.end()){
            mapIter->second++;
        }else {
            classMap.insert(std::pair<int,float>(classe,1));
        }
    }
    
    //        for(std::map<int, float>::iterator it=classMap.begin(); it!=classMap.end(); ++it){
    //            it->second/=indexes.size();
    //        }
    
    
}


//Calcola la media delle features dato un vettore di indici e la inserisce in una mat passata per riferimento

template <typename _T>
 void RandomTree<_T>::classMean(std::vector<unsigned int> &recordIndexes, TrainingData<_T> &dataSet, 
                               cv::Mat_<double> &mean) {
    mean= cv::Mat_<double>::zeros(dataSet.getFeatureSpaceDimension(),1);
    for(unsigned int dim=0;dim<dataSet.getFeatureSpaceDimension();++dim){
        double sum=0.0;
        for(std::vector<unsigned int>::iterator it=recordIndexes.begin();it!=recordIndexes.end();++it){
            sum+=(double)dataSet.getRecord(*it)(dim);
        }
        mean(dim)=((double)sum/recordIndexes.size());

    }
    
}


//Calcola il coefficiente di rayleigh e le medie delle classi e il vettore di proiezione
template <typename _T>
 void RandomTree<_T>::rayleigh(std::vector<unsigned int> &classRecordIndexes, 
                              std::vector<unsigned int> &otherClassesRecordIndexes, 
                              TrainingData<_T> &dataSet, 
                              cv::Mat_<double> &mean1, cv::Mat_<double> &mean2, 
                              cv::Mat_<double> &w, double &rayleigh) {
    
    unsigned int dimensions=dataSet.getFeatureSpaceDimension();
    
    //Matrici in cui vado a salvare i vari record
    cv::Mat_<_T> x1(dimensions,classRecordIndexes.size());
    cv::Mat_<_T> x2(dimensions,otherClassesRecordIndexes.size());
     cv::Mat_<double> sb(dimensions,dimensions);
     cv::Mat_<double> sws1(dimensions,dimensions);
     cv::Mat_<double> sws2(dimensions,dimensions);
     cv::Mat_<double> sw(dimensions,dimensions),swi(dimensions,dimensions);
     
     cv::Mat_<double> ray;
     cv::Mat_<double> x(dimensions,1);
    unsigned int col=0;
    for (std::vector<unsigned int>::iterator it=classRecordIndexes.begin(); it!=classRecordIndexes.end(); ++it) {
        dataSet.getRecordMat(*it,x);
        x1.col(col)=x;
//        for(unsigned int dim=0;dim<dimensions;++dim){
//            x1(dim,col)=dataSet.getRecord(*it)(dim);
//            
//        }
        col++;
    }
    col=0;
    for (std::vector<unsigned int>::iterator it=otherClassesRecordIndexes.begin(); it!=otherClassesRecordIndexes.end(); ++it) {
        dataSet.getRecordMat(*it,x);
        x2.col(col)=x;
//        for(unsigned int dim=0;dim<dataSet.getFeatureSpaceDimension();++dim){
//            x2(dim,col)=dataSet.getRecord(*it)(dim);
//        }
        col++;
    }

    classMean(classRecordIndexes, dataSet, mean1);
    classMean(otherClassesRecordIndexes,dataSet,mean2);

    
    //double rayleigh;
    
    cv::mulTransposed(mean2,sb, false,mean1,1,-1); //calcolo Sb
    cv::mulTransposed(x1, sws1, false,mean1,1,-1); //calcolo s1^
    cv::mulTransposed(x2, sws2, false,mean2,1,-1); //calcolo s2^
    sw=sws1+sws2; //Calcolo Sw
     //invert(sw, swi,CV_SVD); //PseudoInversa di Sw
     
    w=sw.inv(CV_SVD)*(mean2-mean1); //Vettore su cui proietto i dati
     //w=((sw.t()*sw).inv()*sw.t())*(mean2-mean1);
     //w=w/cv::norm(w, cv::NORM_L2); //Versore su cui proietto i dati //spesso in underflow
//     double norm=cv::norm(w,cv::NORM_L2);
//     if (norm!=0) {
//         w=w/norm;
//     }
     ray=(w.t()*sb*w)/(w.t()*sw*w);

    //Stampo il vettore di proiezione per controllo
    //    int c=0;
    //    
    //    cout<<w.rows<<"x"<<w.cols<<endl;
    //    for(Mat_<double>::iterator it=w.begin();it!=w.end();++it){
    //        cout<<*it<<"\t";
    //        c++;
    //        if(c==w.cols){
    //            cout<<endl;
    //            c=0;
    //        }
    //    }
    rayleigh=ray(0);
     //std::cout<<"Rayleigh ratio "<<rayleigh<<std::endl;
    
}

//Calcola un numero casuale tra due valori
template <typename _T>  _T RandomTree<_T>::thresholdRandom(_T valMin, _T valMax)
{
	_T threshold;
	if(valMin == valMax)
		threshold = valMax;
	else
		threshold = (_T) (valMin + rand() % (int)(valMax-valMin));
	
	//Controllo 
    
    //std::cout<<"valMin: "<<valMin<<" valMax: "<<valMax;
    //std::cout << "La soglia e: " << threshold <<std::endl;
    
	//Fine controllo
    return threshold;
}


template <typename _T> 
 void RandomTree<_T>::splitFisher(TreeNode<_T> *node, TrainingData<_T> &dataSet){
    
    /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
     senso proseguire poiche' quello stesso record sarebbe riportato 
     in cascata in un nodo figlio */
    if(node->getTrainingRecordIndexesSize() == 1) {
        node->setIsLeaf(true);
        return;
    }
    /** Se il nodo genitore ha 1 sola classe al suo interno non ha
     senso proseguire */
    std::map<int,float> cMap;
    getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
    if(cMap.size()<=1){
        node->setIsLeaf(true);
        return;
    }
    
    //Ottengo le classi del nodo
    std::map<int, float> classMap;
    getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,classMap);
     
    double r=0,currentr=0;
    unsigned int dimensions=dataSet.getFeatureSpaceDimension();
    cv::Mat_<double> currentmean1=cv::Mat_<double>::zeros(dimensions,1),currentmean2=cv::Mat_<double>::zeros(dimensions,1);
    cv::Mat_<double> sumOfMeans=cv::Mat_<double>::zeros(dimensions,1),meansAverage=cv::Mat_<double>::zeros(dimensions,1);;
    cv::Mat_<double> w(dimensions,1), currentw(dimensions,1);
    std::vector<unsigned int> mainClassIndexes, otherClassesIndexes;
    //Confronto ogni classe con tutte le altre per massimizzare rayleigh
    for(std::map<int,float>::iterator it=classMap.begin();it!=classMap.end();++it){
        mainClassIndexes.clear();
        otherClassesIndexes.clear();
       
        node->getClassFilteredRecords(it->first, dataSet, mainClassIndexes,otherClassesIndexes);
        rayleigh(mainClassIndexes, otherClassesIndexes, dataSet, currentmean1, currentmean2, currentw, currentr);

        if(currentr!=currentr || currentr==0)
            continue;
        if (currentr>r) {
            r=currentr;
            currentw.copyTo(w);
            sumOfMeans=currentmean1+currentmean2;
        }
        
    }
    meansAverage=sumOfMeans/2;
    //Imposto il vettore di proiezione per il nodo
    node->setDirection(w);
    //Imposto la media delle medie per il nodo, serve per la classificazione
    node->setMeansAverage(meansAverage);

        
    /* Allocating left and right child */
    node->setLeftChild(new TreeNode<_T>);
    node->setRightChild(new TreeNode<_T>);
    
    /* Filtro i record del nodo root per il nodo di sinistra */
    node->getFisherFilteredRecords(
                                   w,
                                   meansAverage,
                                   0,
                                   dataSet,
                                   node->leftChild()->getTrainingRecordIndexes());
    
    
    node->getFisherFilteredRecords(
                                   w,
                                   meansAverage,
                                   1,
                                   dataSet,
                                   node->rightChild()->getTrainingRecordIndexes());
    
    
    if(node->leftChild()->getTrainingRecordIndexesSize() == 0 
       || node->rightChild()->getTrainingRecordIndexesSize()==0){
        delete  node->leftChild();
        delete  node->rightChild();
        node->setIsLeaf(true);
        return;
    }
    
    /* Non puÚ accadere che ENTRAMBI i figli siano nulli (a meno che non sia una foglia), 
     * poichË il padre non puÚ essere vuoto e un split binario deve necessariamente spostare 
     * i dati o a destra o a sinistra*/
    
    if(node->leftChild() != 0)
        this->splitFisher(node->leftChild(), dataSet);
    
    if(node->rightChild() != 0)
        this->splitFisher(node->rightChild(), dataSet);
    
}
template <typename _T> 
 void RandomTree<_T>::splitFisher(TreeNode<_T> *node, TrainingData<_T> &dataSet, int height){
    if(height==0){
        node->setIsLeaf(true);
        return;
    }
    /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
     senso proseguire poiche' quello stesso record sarebbe riportato 
     in cascata in un nodo figlio */
    if(node->getTrainingRecordIndexesSize() == 1) {
        node->setIsLeaf(true);
        return;
    }
    std::map<int,float> cMap;
    getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
    if(cMap.size()<=1){
        node->setIsLeaf(true);
        return;
    }
    std::map<int, float> classMap;
    getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,classMap);
    
    double r=0,currentr=0;
    unsigned int dimensions=dataSet.getFeatureSpaceDimension();
    cv::Mat_<double> currentmean1=cv::Mat_<double>::zeros(dimensions, 1),
                     currentmean2=cv::Mat_<double>::zeros(dimensions, 1);
    cv::Mat_<double> sumOfMeans=cv::Mat_<double>::zeros(dimensions, 1),
                     meansAverage=cv::Mat_<double>::zeros(dimensions,1);;
    cv::Mat_<double> w(dimensions,1), currentw(dimensions,1);
    
    //Confronto ogni classe con tutte le altre per massimizzare rayleigh
    for(std::map<int,float>::iterator it=classMap.begin();it!=classMap.end();++it){
        std::vector<unsigned int> mainClassIndexes, otherClassesIndexes;
        
        node->getClassFilteredRecords(it->first, dataSet, mainClassIndexes,otherClassesIndexes);
        rayleigh(mainClassIndexes, otherClassesIndexes, dataSet, currentmean1, currentmean2, currentw, currentr);
        
        if (currentr>r) {
            r=currentr;
            currentw.copyTo(w);
            sumOfMeans=currentmean1+currentmean2;
        }
        
    }
    meansAverage=sumOfMeans/2;
    //Imposto il vettore di proiezione per il nodo
    node->setDirection(w);
    //Imposto la media delle medie per il nodo, serve per la classificazione
    node->setMeansAverage(meansAverage);
    
    
    /* Allocating left and right child */
    node->setLeftChild(new TreeNode<_T>);
    node->setRightChild(new TreeNode<_T>);
    
    /* Filtro i record del nodo root per il nodo di sinistra */
    node->getFisherFilteredRecords(
                                   w,
                                   meansAverage,
                                   0,
                                   dataSet,
                                   node->leftChild()->getTrainingRecordIndexes());
    
    
    node->getFisherFilteredRecords(
                                   w,
                                   meansAverage,
                                   1,
                                   dataSet,
                                   node->rightChild()->getTrainingRecordIndexes());
    
    
    if(node->leftChild()->getTrainingRecordIndexesSize() == 0 
       || node->rightChild()->getTrainingRecordIndexesSize()==0){
        delete  node->leftChild();
        delete  node->rightChild();
        node->setIsLeaf(true);
        return;
    }
    

    
    /* Non puÚ accadere che ENTRAMBI i figli siano nulli (a meno che non sia una foglia), 
     * poichË il padre non puÚ essere vuoto e un split binario deve necessariamente spostare 
     * i dati o a destra o a sinistra*/
    
    if(node->leftChild() != 0)
        this->splitFisher(node->leftChild(), dataSet, height-1);
    
    if(node->rightChild() != 0)
        this->splitFisher(node->rightChild(), dataSet, height-1);
    
}


////Calcolo la soglia che minimizza la differenza del numero di record che andranno nel figio di sinistra e quello di destra.
////Restituisco un pair <distanza minima, valore di soglia>
////
//template <typename _T>
//void RandomTree<_T>::minDiffSplitValue(TreeNode<_T> *node, TrainingData<_T> &dataSet, unsigned int splitIndex, _T valMin, _T valMax, std::pair<float, _T> &minDistElem, float &avgDist) {
//    //Calcolo la distanza tra il numero di features che andranno a destra e sinistra.
//    std::vector<std::pair<unsigned int,_T> > distances;
//
//    std::vector<_T> rvals(getFeatureSpaceDimension());
//    unsigned int totchildrecords=-1;
//    std::vector<unsigned int> recordIndexes=node->getTrainingRecordIndexes();
//    std::map<int, float> classMap;
//    getNodeClassMap(recordIndexes,dataSet,classMap);
//    int maxClass=max_element(classMap.begin(),classMap.end())->first;
//    int maxClassNum=max_element(classMap.begin(),classMap.end(),classProbComparator)->second;
//    avgDist=0;
//    float commonClasses=0;
//    //std::cout<<"Classe principale "<<maxClass<<std::endl;
//    for(std::vector<unsigned int>::iterator it=recordIndexes.begin();it!=recordIndexes.end();++it){
//        for (unsigned int j=0; j<getFeatureSpaceDimension(); ++j) {
//            rvals.push_back(dataSet.getRecord(*it).getFeature(j));
//        }
//    }
//    sort(rvals.begin(),rvals.end());
//    //reverse(rvals.begin(),rvals.end());
//    //    _T i=rvals[0];
//    //    for(int k=1; k<rvals.size();++k){
//    //        if(rvals[k] == rvals[k-1])
//    //            continue;
//    //        i=rvals[k];
//    //        std::cout<<"Soglia "<<i<<std::endl;
//    
//    for(_T i=valMin; i<=valMax;++i){
//        std::vector<unsigned int> lindexes, rindexes;
//        lindexes.clear();
//        rindexes.clear();
//        node->getFilteredRecords(splitIndex, i,0,dataSet,lindexes);
//        node->getFilteredRecords(splitIndex, i,1,dataSet,rindexes);
//        totchildrecords=node->getTrainingRecordIndexesSize();
//        std::map<int, float> lClassMap, rClassMap;
//        getNodeClassMap(lindexes,dataSet,lClassMap);
//        getNodeClassMap(rindexes,dataSet,rClassMap);
//        if (lClassMap.size()==0 || rClassMap.size()==0 
//            || lClassMap.size()/rClassMap.size()!=1
//            //||abs(lClassMap.size()-rClassMap.size())>=classMap.size()/2
//            //||abs(lindexes.size()-rindexes.size())>=totchildrecords/2
//            ) {//non farei lo split
//            continue;
//        }
//        float nPrincL=0, nPrincR=0;
//        if(lClassMap.find(maxClass)!=lClassMap.end())
//            nPrincL=lClassMap.find(maxClass)->second;
//        //nPrincL/=maxClassNum;
//        if(rClassMap.find(maxClass)!=rClassMap.end())
//            nPrincR=rClassMap.find(maxClass)->second;
//        //nPrincR/=maxClassNum;
//        //if(std::abs((int)nPrincL-(int)nPrincR)>totchildrecords/4)
//        // continue;
//        //std::cout<<"Il nodo ha "<<classMap.size()<<" classi della classe "<<maxClass<<" a sinistra vanno "<<nPrincL<<" elementi della classe principale e ci saranno "<<lClassMap.size()<<" classi"<<"a destra "<< nPrincR<<" elementi della classe principale e ci saranno "<<rClassMap.size()<<" classi"<<std::endl;
//        //std::cout<<"Node classes: "<<classMap.size()<<" left classes: "<<lClassMap.size()<<" right classes "<<rClassMap.size()<<std::endl;
//        
//        
//        commonClasses=0;
//        for(std::map<int, float>::iterator liter=lClassMap.begin();liter!=lClassMap.end();++liter){
//            std::map<int, float>::iterator riter=rClassMap.find(liter->first);
//            if (riter!=rClassMap.end()) {
//                commonClasses++;
//            }
//        }
//        commonClasses/=(classMap.size());
//        //        if(std::abs((int)(nPrincL-nPrincR))>totchildrecords/2 && totchildrecords>2)
//        //            continue;
//        //std::cout<<"Common classes: "<<commonClasses<<std::endl;
//        //float distance=commonClasses/classMap.size();
//        //unsigned int distance=std::abs((int)(lindexes.size()-rindexes.size()));
//        float distance=0;
//        //if(commonClasses!=0){//&& std::abs((int)(nPrincL-nPrincR))!=0){
//        distance=std::abs((int)(nPrincL-nPrincR))/commonClasses;
//        //}
//        //        else if (commonClasses!=0 && std::abs((int)(nPrincL-nPrincR))==0){
//        //            distance=1/commonClasses;
//        //        }
//        //        else if (commonClasses==0 && std::abs((int)(nPrincL-nPrincR))==0){
//        //            distance=totchildrecords-std::abs((int)(nPrincL-nPrincR));
//        //        }
//        //unsigned int distance=std::abs((int)(lClassMap.size()-rClassMap.size()));
//        distances.push_back(std::pair<float , _T>(distance,i));
//        //avgDist+=distance;
//    }
//    //std::cout<<std::endl;
//    //sort(distances.begin(),distances.end(),distComparator<_T>); //meglio di no perche' il sort e' lento 
//    //minDistElem=distances.back();
//    if (!distances.empty()) {
//        //sort(distances.begin(),distances.end(),distComparator);
//        minDistElem=*max_element(distances.begin(), distances.end(),distComparator<_T>); 
//        //minDistElem=distances.back();
//        //avgDist/=distances.size();
//        //std::cout<<"Min common classes element "<<minDistElem.second<<std::endl;
//    }else {
//        return;
//    }
//    avgDist=maxClassNum;
//    //if(commonClasses!=0)
//    //avgDist/=commonClasses;
//    distances.clear();
//    rvals.clear();
//    
//}

template <typename _T>
 void RandomTree<_T>::doSplit(TreeNode<_T> *node, TrainingData<_T> &dataSet, unsigned int method) {
    unsigned int spaceDim=dataSet.getFeatureSpaceDimension();
    //Comincio a ricercare il miglior split
    bool badSplit = true;
    unsigned int countRep=0;
    while(badSplit && countRep<=spaceDim) {
        //Calcolo l'indice della feature da usare con lo split in modo random
        unsigned int splitIndex=(rand()%spaceDim);
        
        DataSetIterator it=dataSet.begin();
        _T valMin;
        _T valMax;
        
        
        //Ottengo i valori massimi e minimo delle features del dataset in una data posizione all'interno dei record
        dataSet.getMinMaxValAtFeatureIndex(splitIndex,valMin,valMax);
        
        //valMin = dataSet.getMinValAtFeatureIndex(splitIndex);
        //valMax = dataSet.getMaxValAtFeatureIndex(splitIndex);
        
        _T splitValue=0;
        
        if(valMin==valMax){//Dati allineati
            badSplit=true;
            continue;
            countRep++;
        }
        switch (method) {
            case THRESH_RANDOM:
                //Soglia random
                splitValue=thresholdRandom(valMin,valMax);
                break;
//            case THRESH_MINDIFF:
//                //fisherSplit(node,dataSet,splitIndex);
//                //Calcolo la distanza tra il numero di features che andranno a destra e sinistra.
//                
//                std::pair<float, _T> minDistElem;
//                float avgDist;
//                
//                
//                
//                minDiffSplitValue(node, dataSet, splitIndex, valMin, valMax, minDistElem,avgDist);
//                //std::cout<<"minDistElem.first= "<<minDistElem.first<<" avg "<<avgDist<<std::endl;
//                
//                if((
//                    //                                (minDistElem.first!=0 
//                    //                                 && 
//                    //                minDistElem.first<avgDist/10
//                    //                                 )
//                    //                                || 
//                    (minDistElem.first==0 )
//                    )
//                   && countRep<spaceDim
//                   ){
//                    //                
//                    countRep++;
//                    badSplit=true;
//                    continue;
//                }
//                splitValue=minDistElem.second;
//                //std::cout<<splitValue;
//                
//                break;
        }
        
        
        
        node->setThreshold(splitValue);
        node->setSplitFeatureIndex(splitIndex);
        
        node->setLeftChild(new TreeNode<_T>());
        node->setRightChild(new TreeNode<_T>());
        
        //Imposto i record che andranno nel nodo figlio a sinistra
        node->getFilteredRecords(	splitIndex,
                                 splitValue,
                                 0,
                                 dataSet,
                                 node->leftChild()->getTrainingRecordIndexes());
        
        if(node->leftChild()->getTrainingRecordIndexes().size() == 0) {
            
            /* Elimino i nodi figli*/
            //Il figlio a sinistra non ha record
            delete  node->leftChild();
            delete  node->rightChild();
            
            
            badSplit = true;
            /* Controlla uno split alternativo */
            continue;
        }
        
        //Imposto i record che andranno nel nodo figlio a destra
        node->getFilteredRecords(	splitIndex,
                                 splitValue,
                                 1,
                                 dataSet,
                                 node->rightChild()->getTrainingRecordIndexes());
        
        if(node->rightChild()->getTrainingRecordIndexes().size() == 0) {
            
            /*Elimino i nodi figli*/
            delete  node->leftChild();
            delete  node->rightChild();
            
            
            badSplit = true;
            /* Controlla uno split alternativo */
            continue;
        }
        
        //Qui ho trovato uno split corretto
        badSplit=false;
    }
    
}

template <typename _T>
 void RandomTree<_T>::giniIndexSplit(TreeNode<_T> *node, TrainingData<_T> &dataSet, unsigned int &splitIndex, _T &splitValue , unsigned int method) {
    
    /* Inizializzo il gini score ad un valore alto che sicuramente non sar‡ superato
     durante il calcolo del newGiniScore */
	float currentGiniScore = 9999;
	float newGiniScore;
	
	unsigned int currentBestRecordIndex;
    
	/* Numero di record che dopo lo split selezionato ricadono nel sottoalbero di sinistra */
	unsigned int leftSplitRecords;
	/* Numero di record che dopo lo split selezionato ricadono nel sottoalbero di destra */
	unsigned int rightSplitRecords;
    
	float leftSplitRatio, rightSplitRatio;
	float leftScore, rightScore;
    
	/* Mappa delle occorenze delle classi */
    std::map<unsigned int, unsigned int> classCount;
	
	/* Loop su tutte le dimensioni dello spazio di input. Una possibile variante con il gini index
     (che permetterebbe di usare una foresta di alberi senza che ciascuno di essi dia la stessa risposta
     Ë quella di randomizzare la scelta delle dimensioni per la verifica del giniIndex. Si verificano
     solo un tot di dimensioni scelte in modo random */
    
#ifdef RANDOMIZE_GINI_INDEX
    
	for(unsigned int splitDim = 0; splitDim < this->getFeatureSpaceDimension()/8; splitDim++) {
		unsigned int prevSplitDim = splitDim;
		splitDim = (unsigned int)(rand() % (this->getFeatureSpaceDimension()));
#else
        
        for(unsigned int splitDim = 0; splitDim < this->getFeatureSpaceDimension(); splitDim++) {
            
#endif
            
            
            
            /* Ritorna un vettore con tutti i valori della feature alla dimensione splitDim per ciascun record 
             presente nel nodo genitore  (aggiunge anche il valore della classe per il calcolo del gini index).
             Ordina il vettore in base al valore delle feature */
            
            std::vector<std::pair<_T,unsigned int> > featureValueOnDimensionD;
            
            dataSet.getFeatureValueOnDimensionD(	featureValueOnDimensionD, 
                                                node->getTrainingRecordIndexes(), 
                                                splitDim);
            
            sort(featureValueOnDimensionD.begin(), featureValueOnDimensionD.end(), featureComparator<_T>);
            
            /* Valuto tutti i possibili split su ogni record del nodo, tralasciando i casi in cui due 
             record adiacenti hanno lo stesso valore (non ha senso fare lo split tra questi). Se tutti i
             record hanno lo stesso valore la scelta di quella dimensione per fare lo split non mi porta
             ad alcun risultato e procederÚ oltre nello spazio di input. L'unico problema che si potrebbe 
             presentare Ë se TUTTI i record del nodo rappresentano lo STESSO punto, allora in quel caso
             non si trova alcuno split attuabile. Eventualit‡ remota... (?) */
            
            for(unsigned int recordIndex=1; recordIndex < featureValueOnDimensionD.size(); recordIndex++){
                
                if(featureValueOnDimensionD[recordIndex].first == featureValueOnDimensionD[recordIndex-1].first)
                    continue;
                
                leftSplitRecords = recordIndex;
                rightSplitRecords =  featureValueOnDimensionD.size() - recordIndex;
                
                /* Percentuale dei record che ricadrebbero nei due sottonodi */
                leftSplitRatio  = (float)leftSplitRecords / (float)node->getTrainingRecordIndexesSize();
                rightSplitRatio = (float)rightSplitRecords / (float)node->getTrainingRecordIndexesSize();
                
                /*********************************************************************************************
                 * Valuto per ogni classe il numero di occorrenze nel sottoalbero di sinistra  
                 *********************************************************************************************/
                classCount.clear();
                for(unsigned int recordCounter=0; recordCounter < recordIndex; recordCounter++) {
                    unsigned int c = featureValueOnDimensionD[recordCounter].second;
                    
                    if(classCount.find(c) != classCount.end()) {
                        classCount[c]++;
                    } else {
                        classCount.insert(std::pair<unsigned int, unsigned int>(c,0));
                    }
                }
                
                /* Calcolo per ogni classe esistente lo squared gini ratio */
                
                float temp = 0;
                float squaredGiniRatio;
                
                std::map<unsigned int, unsigned int>::iterator it = classCount.begin();
                
                for(;it!=classCount.end();it++) {
                    squaredGiniRatio = ((float)it->second / (float)leftSplitRecords) * ((float)it->second / (float)leftSplitRecords) ;
                    temp = temp + squaredGiniRatio;
                }
                
                leftScore = 1 - temp;
                leftScore = leftSplitRatio * leftScore;
                
                /*********************************************************************************************
                 * Valuto per ogni classe il numero di occorrenze nel sottoalbero di destra  
                 *********************************************************************************************/
                classCount.clear();
                for(unsigned int recordCounter=recordIndex; recordCounter < featureValueOnDimensionD.size(); recordCounter++) {
                    
                    unsigned int c = featureValueOnDimensionD[recordCounter].second;
                    
                    if(classCount.find(c) != classCount.end()) {
                        classCount[c]++;
                    } else {
                        classCount.insert(std::pair<unsigned int, unsigned int>(c,0));
                    }
                }
                
                /* Per ogni classe calcolo lo squared gini ratio */
                temp = 0;
                squaredGiniRatio = 0;
                
                for(it = classCount.begin();it!=classCount.end();it++) {
                    squaredGiniRatio = ((float)it->second / (float)rightSplitRecords) * ((float)it->second / (float)rightSplitRecords) ;
                    temp = temp + squaredGiniRatio;
                }
                
                rightScore = 1 - temp;
                rightScore = rightSplitRatio * rightScore;
                
                
                
                newGiniScore = leftScore + rightScore;
                if(newGiniScore < currentGiniScore) {
                    currentGiniScore = newGiniScore;
                    splitIndex = splitDim;
                    currentBestRecordIndex = recordIndex;
                    splitValue =  
                    (featureValueOnDimensionD[recordIndex].first - featureValueOnDimensionD[recordIndex-1].first)/2 +
                    featureValueOnDimensionD[recordIndex-1].first;
                }
                
            }
            
            
#ifdef RANDOMIZE_GINI_INDEX
            /* Prima di proseguire nel ciclo riporto il valore di splitDim a quello usato
             all'interno del ciclo */
            splitDim = prevSplitDim;
#endif
            
            
        }
    }
    
    
    
    
    template <typename _T> 
    void RandomTree<_T>::splitGini(TreeNode<_T> *node, TrainingData<_T> &dataSet){
        
        /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
         senso proseguire poiche' quello stesso record sarebbe riportato 
         in cascata in un nodo figlio */
        if(node->getTrainingRecordIndexesSize() == 1) {
            node->setIsLeaf(true);
            return;
        }
        std::map<int,float> cMap;
        getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
        if(cMap.size()<=1){
            node->setIsLeaf(true);
            return;
        }
        
        
        unsigned int splitIndex=0;
        _T splitValue=0;
        
        giniIndexSplit(node, dataSet, splitIndex, splitValue, THRESH_GINI);
        
        node->setSplitFeatureIndex(splitIndex);
        node->setThreshold(splitValue);
        
        /* Allocating left and right child */
        node->setLeftChild(new TreeNode<_T>);
        node->setRightChild(new TreeNode<_T>);
        
        /* Filtro i record del nodo root per il nodo di sinistra */
        node->getFilteredRecords(	splitIndex,
                                 splitValue,
                                 0,
                                 dataSet,
                                 node->leftChild()->getTrainingRecordIndexes());
        
        
        node->getFilteredRecords(	splitIndex,
                                 splitValue,
                                 1,
                                 dataSet,
                                 node->rightChild()->getTrainingRecordIndexes());
        
        
        if(node->leftChild()->getTrainingRecordIndexesSize() == 0 
           || node->rightChild()->getTrainingRecordIndexesSize()==0){
            delete  node->leftChild();
            delete  node->rightChild();
            node->setIsLeaf(true);
            return;
        }
        
        /* Non puÚ accadere che ENTRAMBI i figli siano nulli (a meno che non sia una foglia), 
         * poichË il padre non puÚ essere vuoto e un split binario deve necessariamente spostare 
         * i dati o a destra o a sinistra*/
        
        if(node->leftChild() != 0)
            this->splitGini(node->leftChild(), dataSet);
        
        if(node->rightChild() != 0)
            this->splitGini(node->rightChild(), dataSet);
        
    }
    template <typename _T> 
     void RandomTree<_T>::splitGini(TreeNode<_T> *node, TrainingData<_T> &dataSet, int height){
        if(height==0){
            node->setIsLeaf(true);
            return;
        }
        /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
         senso proseguire poiche' quello stesso record sarebbe riportato 
         in cascata in un nodo figlio */
        if(node->getTrainingRecordIndexesSize() == 1) {
            node->setIsLeaf(true);
            return;
        }
        std::map<int,float> cMap;
        getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
        if(cMap.size()<=1){
            node->setIsLeaf(true);
            return;
        }
        
        unsigned int splitIndex = 0;
        _T splitValue = 0;
        
        giniIndexSplit(node, dataSet, splitIndex, splitValue, THRESH_GINI);
        
        node->setSplitFeatureIndex(splitIndex);
        node->setThreshold(splitValue);
        
        /* Allocating left and right child */
        node->setLeftChild(new TreeNode<_T>);
        node->setRightChild(new TreeNode<_T>);
        
        /* Filtro i record del nodo root per il nodo di sinistra */
        node->getFilteredRecords(	splitIndex,
                                 splitValue,
                                 0,
                                 dataSet,
                                 node->leftChild()->getTrainingRecordIndexes());
        
        
        node->getFilteredRecords(	splitIndex,
                                 splitValue,
                                 1,
                                 dataSet,
                                 node->rightChild()->getTrainingRecordIndexes());
        
        
        if(node->leftChild()->getTrainingRecordIndexesSize() == 0 
           || node->rightChild()->getTrainingRecordIndexesSize()==0){
            delete  node->leftChild();
            delete  node->rightChild();
            node->setIsLeaf(true);
            return;
        }
        
        /* Non puÚ accadere che ENTRAMBI i figli siano nulli (a meno che non sia una foglia), 
         * poichË il padre non puÚ essere vuoto e un split binario deve necessariamente spostare 
         * i dati o a destra o a sinistra*/
        
        if(node->leftChild() != 0)
            this->splitGini(node->leftChild(), dataSet, height-1);
        
        if(node->rightChild() != 0)
            this->splitGini(node->rightChild(), dataSet, height-1);    
    }
//    template <typename _T>
//    void RandomTree<_T>::splitMinDiff(TreeNode<_T> *node, TrainingData<_T> &dataSet) {
//        
//        /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
//         senso proseguire poiche' quello stesso record sarebbe riportato 
//         in cascata in un nodo figlio */
//        if(node->getTrainingRecordIndexesSize() == 1) {
//            node->setIsLeaf(true);
//            return;
//        }
//        std::map<int,float> cMap;
//        getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
//        if(cMap.size()<=1){
//            node->setIsLeaf(true);
//            return;
//        }
//        
//        
//        //Eseguo lo split ottimizzando la distribuzione dei record tra i figli
//        doSplit(node, dataSet, THRESH_MINDIFF);
//        
//        
//        if(node->leftChild()->getTrainingRecordIndexesSize() == 0 
//           || node->rightChild()->getTrainingRecordIndexesSize()==0){
//            node->setIsLeaf(true);
//            return;
//        }
//        
//        if(node->leftChild() != NULL)
//            this->splitMinDiff(node->leftChild(), dataSet);
//        
//        if(node->rightChild() != NULL)
//            this->splitMinDiff(node->rightChild(), dataSet);
//    }
//    
//    
//    template <typename _T>
//    void RandomTree<_T>::splitMinDiff(TreeNode<_T> *node, TrainingData<_T> &dataSet, int height) {
//        if(height==0){
//            node->setIsLeaf(true);
//            return;
//        }
//        
//        /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
//         senso proseguire poiche' quello stesso record sarebbe riportato 
//         in cascata in un nodo figlio */
//        if(node->getTrainingRecordIndexesSize() == 1) {
//            node->setIsLeaf(true);
//            return;
//        }
//        std::map<int,float> cMap;
//        getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
//        if(cMap.size()<=1){
//            node->setIsLeaf(true);
//            return;
//        }
//        
//        doSplit(node, dataSet, THRESH_MINDIFF);
//        
//        
//        if(node->leftChild()->getTrainingRecordIndexesSize() == 0 
//           || node->rightChild()->getTrainingRecordIndexesSize()==0){
//            node->setIsLeaf(true);
//            return;
//        }
//        if(node->leftChild() != NULL)
//            this->splitMinDiff(node->leftChild(), dataSet, height -1);
//        
//        if(node->rightChild() != NULL)
//            this->splitMinDiff(node->rightChild(), dataSet, height-1);
//    }
    
    template <typename _T>
    void RandomTree<_T>::splitRandom(TreeNode<_T> *node, TrainingData<_T> &dataSet) {
        
        /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
         senso proseguire poiche' quello stesso record sarebbe riportato 
         in cascata in un nodo figlio */
        if(node->getTrainingRecordIndexesSize() == 1) {
            node->setIsLeaf(true);
            return;
        }
        std::map<int,float> cMap;
        getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
        if(cMap.size()<=1){
            node->setIsLeaf(true);
            return;
        }
        
        //Eseguo lo split con soglia random
        doSplit(node, dataSet, THRESH_RANDOM);
        
        if(node->leftChild()->getTrainingRecordIndexesSize() == 0 
           || node->rightChild()->getTrainingRecordIndexesSize()==0){
            node->setIsLeaf(true);
            return;
        }
        
        if(node->leftChild() != NULL)
            this->splitRandom(node->leftChild(), dataSet);
        
        if(node->rightChild() != NULL)
            this->splitRandom(node->rightChild(), dataSet);
    }
    
    template <typename _T>
    void RandomTree<_T>::splitRandom(TreeNode<_T> *node, TrainingData<_T> &dataSet, int height) {
        if(height==0){
            node->setIsLeaf(true);
            return;
        }
        
        /** Se il nodo genitore ha 1 solo elemento al suo interno non ha
         senso proseguire poiche' quello stesso record sarebbe riportato 
         in cascata in un nodo figlio */
        if(node->getTrainingRecordIndexesSize() == 1) {
            node->setIsLeaf(true);
            return;
        }
        std::map<int,float> cMap;
        getNodeClassMap(node->getTrainingRecordIndexes(),dataSet,cMap);
        if(cMap.size()<=1){
            node->setIsLeaf(true);
            return;
        }
        
        
        //Eseguo lo split con soglia random
        doSplit(node, dataSet, THRESH_RANDOM);
        
        if(node->leftChild() != NULL)
            this->splitRandom(node->leftChild(), dataSet, height -1);
        
        if(node->rightChild() != NULL)
            this->splitRandom(node->rightChild(), dataSet, height-1);
    }
    
    /**
     * Funzione che implementa l'algoritmo di training del random tree
     */
    
    template <typename _T>  void RandomTree<_T>::trainTree(TrainingData<_T> &dataSet, int height, unsigned int method) {
        
        /* Imposta la dimensione dello spazio di input */ 
        this->setFeatureSpaceDimension(dataSet.getFeatureSpaceDimension());
        this->_DataSet=dataSet;
        //Creo un nuovo nodo
        this->_root = new TreeNode<_T>();
        
        /* Inserisco nel nodo root gli indici di tutti gli elmenti di
         training */
        
        DataSetIterator it=dataSet.begin();
        for(; it != dataSet.end(); it++) {
            _root->pushRecordIndex(it-dataSet.begin());
        }
        
        /* Comincio la fase di split ricorsivo */
        
        switch (method) {
            case THRESH_FISHER:
                //Split con distribuzione ottimale dei record
                //this->splitMinDiff(this->_root, dataSet, height-1);
                this->splitFisher(this->_root,dataSet,height-1);
                break;
            case THRESH_RANDOM:
                //Split con soglia random
                this->splitRandom(this->_root, dataSet, height-1);
                break;
            case THRESH_GINI:
                this->splitGini(this->_root,dataSet, height-1);
                break;
            default:
                //Split con soglia random
                this->splitRandom(this->_root, dataSet, height-1);
                break;
        }
        
        
        
    }
    
    /**
     * Funzione che implementa l'algoritmo di training del random tree
     */
    
    template <typename _T>  void RandomTree<_T>::trainTree(TrainingData<_T> &dataSet, unsigned int method) {
        
        /* Imposta la dimensione dello spazio di input */ 
        this->setFeatureSpaceDimension(dataSet.getFeatureSpaceDimension());
        this->_DataSet=dataSet;
        //Creo un nuovo nodo
        this->_root = new TreeNode<_T>();
        
        /* Inserisco nel nodo root gli indici di tutti gli elmenti di
         training */
        
        DataSetIterator it=dataSet.begin();
        for(; it != dataSet.end(); it++) {
            _root->pushRecordIndex(it-dataSet.begin());
        }
        
        /* Comincio la fase di split ricorsivo */
        
        switch (method) {
            case THRESH_FISHER:
                //Split con distribuzione ottimale dei record
                //this->splitMinDiff(this->_root, dataSet);
                this->splitFisher(this->_root,dataSet);
                break;
            case THRESH_RANDOM:
                //Split con soglia random
                this->splitRandom(this->_root, dataSet);
                break;
                
                //        case THRESH_IG:
                //            this->splitIG(this->_root, dataSet);
                //            break;
                
            case THRESH_GINI:
                this->splitGini(this->_root,dataSet);
                break;
            default:
                //Split con soglia random
                this->splitRandom(this->_root, dataSet);
                break;
        }
        
        
        
    }
    
    
    //Eseguo la classificazione di un record
    template <typename _T>  unsigned int RandomTree<_T>::classify(MyRecord<_T> &record) {
        
        
        TreeNode<_T> *currNode = _root;
        
        
        while(currNode->isLeaf() == false) {
            
            /* Naviga l'albero seguendo i vari split */
            if(record(currNode->getSplitFeatureIndex()) <= currNode->getThreshold()) {
                
                /* Se il nodo in cui devo scendere e' nullo
                 mi fermo al nodo che ho raggiunto fino a questo momento */
                if(currNode->leftChild() != NULL)
                    currNode = currNode->leftChild(); 
                else 
                    break;
            }else { 
                if(currNode->rightChild() != NULL)
                    currNode = currNode->rightChild(); 
                else 
                    break;
                
            }
        }
        
        /* Dopo aver raggiunto una foglia viene fatto l'istogramma delle classi */
        std::vector<unsigned int>& trainingRecords = currNode->getTrainingRecordIndexes();
        std::vector<unsigned int>::iterator it = trainingRecords.begin();
        
        std::map<unsigned int, unsigned int> classFreq;

        for(;it!= trainingRecords.end(); it++) {
            if(classFreq.find(this->_DataSet.getRecord(*it).getClass()) != classFreq.end()) {
                classFreq[this->_DataSet.getRecord(*it).getClass()] ++;
            } else {
                classFreq[this->_DataSet.getRecord(*it).getClass()]= 1;
            }
        }
        unsigned int classe=0, count=0;
        for(std::map<unsigned int, unsigned int>::iterator it=classFreq.begin(); it != classFreq.end(); it++)
        {
            if(count < it->second)
            {
                classe = it->first;
                count = it->second;
            }
            
        }
        return classe;  
    }
    
    //Eseguo la classificazione di un record
    template <typename _T>  unsigned int RandomTree<_T>::classifyFisher(cv::Mat_<double> &record) {
        
        
        TreeNode<_T> *currNode = _root;
        
        
        while(currNode->isLeaf() == false) {
            cv::Mat_<double> res,direction(this->_DataSet.getFeatureSpaceDimension(),1),meansAverage(this->_DataSet.getFeatureSpaceDimension(),1);
   
            currNode->getDirection(direction);
            currNode->getMeansAverage(meansAverage);
            //Classificazione secondo fisher
            res=direction.t()*record-direction.t()*(meansAverage);
            double r=res(0);
            /* Naviga l'albero seguendo i vari split */
            if(r <= 0) {
                
                /* Se il nodo in cui devo scendere e' nullo
                 mi fermo al nodo che ho raggiunto fino a questo momento */
                if(currNode->leftChild() != NULL)
                    currNode = currNode->leftChild(); 
                else 
                    break;
            }else { 
                if(currNode->rightChild() != NULL)
                    currNode = currNode->rightChild(); 
                else 
                    break;
                
            }
        }
        
        /* Dopo aver raggiunto una foglia viene fatto l'istogramma delle classi */
        std::vector<unsigned int>& trainingRecords = currNode->getTrainingRecordIndexes();
        std::vector<unsigned int>::iterator it = trainingRecords.begin();
        
        std::map<unsigned int, unsigned int> classFreq;
        for(;it!= trainingRecords.end(); it++) {
            if(classFreq.find(this->_DataSet.getRecord(*it).getClass()) != classFreq.end()) {
                classFreq[this->_DataSet.getRecord(*it).getClass()] ++;
            } else {
                classFreq[this->_DataSet.getRecord(*it).getClass()]= 1;
            }
        }
        unsigned int classe=0, count=0;
        for(std::map<unsigned int, unsigned int>::iterator it=classFreq.begin(); it != classFreq.end(); it++)
        {
            if(count < it->second)
            {
                classe = it->first;
                count = it->second;
            }
            
        }
        return classe;  
    }

    
    //Imposto la dimensionalita' delle features
    template <typename _T> void RandomTree<_T>::setFeatureSpaceDimension(unsigned int featureSpaceDimension){
        this->_featureSpaceDimension = featureSpaceDimension;
    }
    
    //Restituisco la dimensionalita' delle features
    template <typename _T> unsigned int RandomTree<_T>::getFeatureSpaceDimension(){
        return this->_featureSpaceDimension;
    }
    
    //Istanzio le classi con i possibili tipi di dato da utilizzare
    template class RandomTree<int>;
    template class RandomTree<float>;
    template class RandomTree<double>;
    template class RandomTree<unsigned char>;